#ifndef B_BASIS_FUNCTION_H
#define B_BASIS_FUNCTION_H

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
#include "Polynomials.h"
#include <stdexcept>
//计算基函数，一阶导，二阶导
class B_Basis_function{
    public:
    int degree;
    int index;
    std::vector<double> knots;


    public:
    B_Basis_function(){}
    B_Basis_function(int degree,int index,std::vector<double>knots);
    double value(double x);
    double m(double x_value);
    double M(double x_value);
    double checkIndex(int index);

};

B_Basis_function::B_Basis_function(int degree,int index,std::vector<double>knots):degree(degree),index(index),knots(knots){}
    double B_Basis_function::checkIndex(int index) {
    int n = knots.size(); // 节点数
    if (index < 0) {
        // 如果 index 在范围外，返回左侧的扩展值
        return knots[0] + index;
    } else if (index >= n) {
        // 如果 index 在范围外，返回右侧的扩展值
        return knots[n-1] + (index - (n-1)) ;
    } else {
        // 如果 index 在节点范围内，直接返回节点值
        return knots[index];
    }
}
     /*
double B_Basis_function::value(double x){
    if (degree == 0)
    {
        if (( checkIndex(index-1)<=x) && (x < checkIndex(index)))
        {
            return 1;
        }
        else
        {
            return 0;
        }        
    }
    else
    {
        double coef_1 = (double)(x-checkIndex(index-1))/(double)(checkIndex(index+degree-1)-checkIndex(index-1));
        double coef_2 = (double)(checkIndex(index+degree)-x)/(double)(checkIndex(index+degree)-checkIndex(index));
        B_Basis_function B_1(degree-1, index, knots);
        B_Basis_function B_2(degree-1, index+1, knots);
        return (coef_1*B_1.value(x)+coef_2*B_2.value(x));
    }  
}
*/
double B_Basis_function::value(double x) {
    if (degree == 0) {
        double lower_bound = checkIndex(index - 1);
        double upper_bound = checkIndex(index);
       // std::cout << "Degree 0: x = " << x << ", lower_bound = " << lower_bound << ", upper_bound = " << upper_bound << std::endl;
        if (x>=lower_bound && x < upper_bound) {
            return 1;
        } else {
            return 0;
        }
    } else {
        double denom1 = checkIndex(index + degree - 1) - checkIndex(index - 1);
        double denom2 = checkIndex(index + degree) - checkIndex(index);

       // std::cout << "Degree " << degree << ": x = " << x << ", denom1 = " << denom1 << ", denom2 = " << denom2 << std::endl;

        // 添加边界条件检查，确保分母不为零
        if (denom1 == 0 || denom2 == 0) {
            std::cerr << "Error: Denominator is zero in B_Basis_function::value" << std::endl;
            return 0;
        }

        double coef_1 = (x - checkIndex(index - 1)) / denom1;
        double coef_2 = (checkIndex(index + degree) - x) / denom2;

        //std::cout << "Degree " << degree << ": coef_1 = " << coef_1 << ", coef_2 = " << coef_2 << std::endl;

        B_Basis_function B_1(degree - 1, index, knots);
        B_Basis_function B_2(degree - 1, index + 1, knots);

        double result = (coef_1 * B_1.value(x) + coef_2 * B_2.value(x));
       // std::cout << "Degree " << degree << ": result = " << result << std::endl;

        return result;
    }
}
double B_Basis_function::m(double x_value){
    if(degree==0){
        return 0;
    }else{
        
        B_Basis_function b1(degree-1,index,knots);
        B_Basis_function b2(degree-1,index+1,knots);
        double denom1 = checkIndex(index+degree-1) - checkIndex(index-1);
        double denom2 = checkIndex(index+degree) - checkIndex(index);
        if (denom1 != 0 && denom2 != 0) {
            //std::cout<<"一阶导数："<<(degree)*b1.value(x_value)/denom1 - (degree)*b2.value(x_value)/denom2<<endl;
            return (degree)*b1.value(x_value)/denom1 - (degree)*b2.value(x_value)/denom2;
        } else {
            return 0.0;
        }
    }
}
double B_Basis_function::M(double x_value){
    if(degree==0||degree==1){
        return 0;
    }else{
        B_Basis_function b1(degree-2,index,knots);
        B_Basis_function b2(degree-2,index+1,knots);
        double denom1 = (checkIndex(index+degree-1) - checkIndex(index-1))*(checkIndex(index+degree-2)-checkIndex(index-1));
        double denom2 = (checkIndex(index+degree) - checkIndex(index))*(checkIndex(index+degree-1)-checkIndex(index));
        if (denom1 != 0 && denom2 != 0) {
//std::cout<<"二阶导数："<<(degree*(degree-1))*b1.value(x_value)/denom1 - (degree*(degree-1))*b2.value(x_value)/denom2<<endl;
            return (degree*(degree-1))*b1.value(x_value)/denom1 - (degree*(degree-1))*b2.value(x_value)/denom2;
        } else {
            return 0.0;
        }
    }
}
#endif